
stub<-"h"
nimp<-100
alpha<-0.95
jitter<-0.2

print("confirm version")

impdata<-read.csv("baummerge2.csv")


#subcty<-unique(impdata$countrycode)

#subcty<-c("CUB","MEX","CRI","CHL","GBR","CIV","CMR","CAF","TZA","SAU","OMN")

#subcty<-c("CUB","MEX","CIV")

#subcty<-c("MEX","CHL","FRA","CPV","CIV","PRK")

#subcty<-c("USA","CAN","DOM","JAM","CUB","MEX","SLV","CRI","PER","BRA","CHL","GBR","FRA","CHE","ITA","SWE","NOR","CIV","CMR","CAF","TZA","SAU","OMN","MYS","SGP","AUS","NZL")

#subcty<-c("MEX","CIV","CMR")

subcty<-c("ECU","PER","IRL","GRC","SDN","BHR","MNG","NPL","MYS")

subcty<-c("ECU","GRC","MNG","NPL")


viewedlimits<-matrix(c(-5,70,10,110,10,110,-5,50),2,4)
viewedlimits<-t(viewedlimits)


ci.upper<-function(data){
  alpha<-0.95   #local alpha, not global alpha, makes apply easier, and yet cumbersome
  a<-sort(data)
  b<-a[round(length(a)*alpha)] 
#  b<-a[length(a)] 
  return(b)
}

ci.lower<-function(data){
  alpha<-0.95
  a<-sort(data)
  b<-a[round(length(a)*(1-alpha))] 
#  b<-a[1] 
  return(b)
}



period<-c(1959,2000)


yearseq<-(period[1]+1):(period[2]-1)
xyears<-rbind(yearseq,yearseq,NA)
delta<-0.1
xyearsdelta<-rbind(yearseq-delta,yearseq+delta,NA)

l<-1
w<-4
lw<-l*w
par(mfrow=c(l,w)) 




for(i in 1:length(subcty)){

  #allmain<-countrynames[i]

  flag<-impdata$countrycode==subcty[i] & impdata$year>period[1] & impdata$year<period[2]
  tempdata<-impdata[flag, ]
  allmain<-tempdata[1,"countryname"]

#  lowerb<-min(tempdata$femschool,na.rm=TRUE) - 0.5*abs(min(tempdata$femschool,na.rm=TRUE))
#  upperb<-max(tempdata$femschool,na.rm=TRUE)*1.3
#  allylim<-c(lowerb,upperb)

  allylim<-viewedlimits[i,]

 
  r1<-is.na(tempdata$femschool)
  impeddata<-matrix(0,period[2]-period[1]-1,nimp)
  stackdata<-matrix(0,period[2]-period[1]-1,nimp)

  for(j in 1:nimp){
    temp.imp<-read.csv(file=paste(stub,j,".csv",sep=""))
    impflag<-temp.imp$countrycode==subcty[i] & temp.imp$year>period[1] & temp.imp$year<period[2]
    impeddata[,j]<-temp.imp$femschool[impflag]
  }

  for(i in 1:nrow(stackdata)){
    stackdata[i,]<-sort(impeddata[i,])
  }

#  s.upper <-apply(stackdata,1,ci.upper)
#  s.median<-apply(stackdata,1,median)
#  s.lower <-apply(stackdata,1,ci.lower)

   cilr<-c(round(nimp*(1-alpha)),round(nimp*alpha))

   missing<-is.na(tempdata$femschool)
   years.missing<-xyears[1,missing]


   s.upper <-  stackdata[,cilr[1]]
   s.median <- stackdata[,round(nimp/2)]
   s.lower <-  stackdata[,cilr[2]]

#print(cbind(stackdata,s.upper,s.lower))
#stop()


#  # Imputed Values
#  plot((years.missing*matrix(1,length(years.missing),cilr[1]-1))+runif((length(years.missing))*(cilr[1]-1),-jitter,jitter),stackdata[missing,1:(cilr[1]-1)],main=allmain,xlab="Year",ylab="Female Secondary Schooling",ylim=allylim,xlim=period,type="p",col="salmon")
#  par(new=TRUE)
#  plot((years.missing*matrix(1,length(years.missing),nimp-cilr[2]))+runif((length(years.missing))*(nimp-cilr[2]),-jitter,jitter),stackdata[missing,(cilr[2]+1):nimp],main=allmain,xlab="Year",ylab="Female Secondary Schooling",ylim=allylim,xlim=period,type="p",col="salmon")
#  par(new=TRUE)
#  plot((years.missing*matrix(1,length(years.missing),cilr[2]-cilr[1]+1))+runif((length(years.missing))*(cilr[2]-cilr[1]+1),-jitter,jitter),stackdata[missing,cilr[1]:cilr[2]],main=allmain,xlab="Year",ylab="Female Secondary Schooling",ylim=allylim,xlim=period,type="p",col="grey")
#  par(new=TRUE)


  # Confidence Intervals  
  plot(xyears[,r1],rbind(s.upper[r1],s.lower[r1],NA),main=allmain,xlab="Year",ylab="Female Secondary Schooling",ylim=allylim,xlim=period,type="l",col="slategrey")
  par(new=TRUE)
  plot(xyearsdelta[,r1],rbind(s.median[r1],s.median[r1],NA),main=allmain,xlab="Year",ylab="Female Secondary Schooling",ylim=allylim,xlim=period,type="l",col="slategrey")
  par(new=TRUE)

  # Some (m) Imputed Values
  m<-5
  plot((years.missing*matrix(1,length(years.missing),m)),impeddata[missing,31:(30+m)],main=allmain,xlab="Year",ylab="Female Secondary Schooling",ylim=allylim,xlim=period,type="p",col="blue")
  par(new=TRUE)

  # Observed Data
  plot(tempdata$year,tempdata$femschool,main=allmain,xlab="Year",ylab="Female Secondary Schooling",ylim=allylim,xlim=period,pch=16)

  
#  if(round(i/lw)*lw==i){
#    par(ask=TRUE)
#  }else{
#    par(ask=FALSE)
#  }

}

  
